Bloque 4

¿Qué haremos en este bloque?

Cuando trabajamos con diseños complejos no podemos usar métodos estadísticos tradicionales sin ajustar la varianza y los ponderadores muestrales. En este bloque exploraremos modelos de regresión lineal, logística y técnicas avanzadas para análisis de encuestas.

Introducción a los modelos con encuestas

Preparamos los datos

Cargamos las librerías de trabajo:

Click para ver el código

# Para instalar manualmente 
install.packages("tidyverse")
install.packages("rio")
install.packages("haven")
install.packages("survey")
install.packages("srvyr")
install.packages("vtable")
install.packages("tidymodels")
install.packages("texreg")
install.packages("factoextra")
install.packages("NbClust")

# Cargamos las librerías de trabajo 
library("tidyverse")
library("rio")
library("haven")
library("survey")
library("srvyr")
library("vtable")
library("tidymodels")
library("texreg")
library("factoextra")
library("NbClust")

Abrimos los datos:

Click para ver el código

# Hay varios caminos para abrir los datos dependiendo el formato:
data <- rio::import("20240516_enssex_data.rdata") # R

Preparamos las variables:

Click para ver el código

variables <- c("folio_encuesta", "macro_region", "region", 
               "varstrat", # Pseuestrato
               "varunit", # Pseudoconglomerado
               "w_personas_cal", "w_personas_cal_nr_mod_hogar", # Ponderadores
               "fecha", "edad_primera_rsex", "relaciones_sexuales", 
               "p1", "genero_separada", "orientacion_sexual_separado",
               "edad", "fecha_nacimiento_ano", 
               "p34", # Conversaciones sobre sexualidad
               "i_6_p9", # Bienestar con la vida sexual
               "educ", 
               "p247", # Victima de abuso sexual
               "evaluacion_educacion_sexual", 
               "lugar_primera_rsex", 
               "vinculo_primera_rsex")

data <- data |> 
  dplyr::select(all_of(variables)) |> 
  dplyr::mutate(fecha_encuesta=dmy(fecha),
                agno_encuesta=year(fecha_encuesta), 
                edad_val=agno_encuesta-edad, 
                validacion_fecha_nac=if_else(edad_val==fecha_nacimiento_ano, 1, 0)) |> 
  mutate(
    macro_region=factor(macro_region, labels=c("Macrozona Norte",
                                               "Macrozona Centro", 
                                               "Macrozona Sur")),
    region=factor(region, labels=c("Tarapacá", 
                                   "Antofagasta", 
                                   "Atacama", 
                                   "Coquimbo",
                                   "Valparaíso", 
                                   "O'Higgins", 
                                   "Maule", 
                                   "Biobío", 
                                   "La Araucanía", 
                                   "Los Lagos", 
                                   "Aysén", 
                                   "Magallanes", 
                                   "Metropolitana", 
                                   "Los Ríos", 
                                   "Arica y Parinacota", 
                                   "Ñuble")),
    sex=factor(p1, labels=c("Hombre", "Mujer")), 
    genero=factor(genero_separada, labels=c("Masculino",
                                            "Transgénero", 
                                            "No binario", 
                                            "Otro",
                                            "Femenino",
                                            "No responder")),
    vinculo_primera_rsex=factor(vinculo_primera_rsex, labels = c("Pareja", 
                                                       "Amigo/a",
                                                       "Amigo/a con beneficios", 
                                                       "Otro",
                                                       "No responde")),
     lugar_primera_rsex=factor(lugar_primera_rsex, labels = c("Casa propia", 
                                                       "Casa de los padres",
                                                       "Casa de amigos/as", 
                                                       "Hotel",
                                                       "Aire libre",
                                                       "Otro",
                                                       "No responde")),
    cohorte=case_when(
      fecha_nacimiento_ano<1950 ~ 1, 
      fecha_nacimiento_ano>=1950 & fecha_nacimiento_ano<1960 ~ 2,
      fecha_nacimiento_ano>=1960 & fecha_nacimiento_ano<1970 ~ 3,
      fecha_nacimiento_ano>=1970 & fecha_nacimiento_ano<1980 ~ 4,
      fecha_nacimiento_ano>=1980 & fecha_nacimiento_ano<1990 ~ 5,
      fecha_nacimiento_ano>=1990 & fecha_nacimiento_ano<2000 ~ 6,
      fecha_nacimiento_ano>=2000 ~ 7,
      TRUE ~ NA
    ), 
    cohorte=factor(cohorte, labels = c("Cohorte < 1950", 
                                     "Cohorte 1951-1960",
                                     "Cohorte 1961-1970",
                                     "Cohorte 1971-1980",
                                     "Cohorte 1981-1990",
                                     "Cohorte 1991-2000",
                                     "Cohorte 2001-2005")),
    p34=if_else(p34 %in% c(8,9), 3, p34-1),
    p34=if_else(p34==2, 1, p34),
    sexual_convers=factor(p34, labels = c("No", "Si", "No responde")), 
    educ = factor(educ, labels=c("Primaria",
                                 "Secundaria incompleta",
                                 "Secondarian completa",
                                 "Universitario")),
    sexual_abuso = factor(p247, labels=c("Si", "No", "No responde")),
    sexual_educ=factor(evaluacion_educacion_sexual, labels=c("Mal", 
                                                             "Regular", 
                                                             "Bien", 
                                                             "Sin respuesta"
                                                                  )),

    bienestar=if_else(i_6_p9==9, NA, i_6_p9),
    bienestar=as.numeric(bienestar)
  ) |> 
  dplyr::select(varstrat, varunit, w_personas_cal, 
         macro_region, region,
         bienestar,
         edad_primera_rsex, 
         vinculo_primera_rsex, lugar_primera_rsex,
         sex, orientacion_sexual_separado, 
         edad, cohorte, educ,
         sexual_convers,
         sexual_educ, 
         sexual_abuso
         ) |> 
  drop_na()

Revisamos nuestra tabla descriptiva:

Click para ver el código

# Tabla de descriptivos sin ponderar 
st(dplyr::select(data, !c("varstrat", "varunit", "w_personas_cal", "region")),
   #labels=c("Identidad de género", "Grupo Etario", "Nivel educativo", "Víctima de violencia"),
   digits = 3, 
   out="kable", 
   add.median = TRUE,
   fixed.digits = TRUE, 
   simple.kable = FALSE,
   title="",
   numformat = NA) %>%
  kableExtra::kable_styling(latex_options = c("scale_down", "hold_position"),
                            full_width = FALSE, fixed_thead = T)

Variable	N	Mean	Std. Dev.	Min	Pctl. 25	Pctl. 50	Pctl. 75	Max
macro_region	16904
... Macrozona Norte	3261	19.3%
... Macrozona Centro	8457	50.0%
... Macrozona Sur	5186	30.7%
bienestar	16904	4.964	1.997	1	4	6	7	7
edad_primera_rsex	16904	17.739	3.484	10	16	17	19	76
vinculo_primera_rsex	16904
... Pareja	10319	61.0%
... Amigo/a	2164	12.8%
... Amigo/a con beneficios	2267	13.4%
... Otro	2082	12.3%
... No responde	72	0.4%
lugar_primera_rsex	16904
... Casa propia	4812	28.5%
... Casa de los padres	5434	32.1%
... Casa de amigos/as	1686	10.0%
... Hotel	893	5.3%
... Aire libre	1893	11.2%
... Otro	1824	10.8%
... No responde	362	2.1%
sex	16904
... Hombre	5805	34.3%
... Mujer	11099	65.7%
edad	16904	45.104	17.123	18	30	43	58	100
cohorte	16904
... Cohorte < 1950	1172	6.9%
... Cohorte 1951-1960	1986	11.7%
... Cohorte 1961-1970	2693	15.9%
... Cohorte 1971-1980	2806	16.6%
... Cohorte 1981-1990	3139	18.6%
... Cohorte 1991-2000	3814	22.6%
... Cohorte 2001-2005	1294	7.7%
educ	16904
... Primaria	3033	17.9%
... Secundaria incompleta	2365	14.0%
... Secondarian completa	5765	34.1%
... Universitario	5741	34.0%
sexual_convers	16904
... No	11893	70.4%
... Si	4841	28.6%
... No responde	170	1.0%
sexual_educ	16904
... Mal	8167	48.3%
... Regular	4311	25.5%
... Bien	3571	21.1%
... Sin respuesta	855	5.1%
sexual_abuso	16904
... Si	2106	12.5%
... No	13647	80.7%
... No responde	1151	6.8%

Regresión lineal

La regresión lineal ponderada estima la relación entre una variable dependiente continua \(Y\) y un conjunto de predictores \(X_1, X_2, \dots, X_p\):

\[Y_i = \beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2} + \dots + \beta_p X_{ip} + \varepsilon_i\]

Donde:
- \(Y_i\) es la variable dependiente (Ej: edad de inicio sexual).
- \(X_{i1}, X_{i2}, \dots, X_{ip}\) son las variables explicativas (Ej: sexo, nivel educativo, región).
- \(\beta_0\) es la intersección (constante).
- \(\beta_p\) son los coeficientes de regresión.
- \(\varepsilon_i\) es el error aleatorio.

Click para ver el código

# Sin diseño muestral
m1 <- lm(bienestar ~ edad_primera_rsex + sex + orientacion_sexual_separado + 
        edad + cohorte + educ + sexual_convers + sexual_educ, data=data)

# Podemos ver los resultados de distintas maneras:
summary(m1)


Call:
lm(formula = bienestar ~ edad_primera_rsex + sex + orientacion_sexual_separado + 
    edad + cohorte + educ + sexual_convers + sexual_educ, data = data)

Residuals:
   Min     1Q Median     3Q    Max 
-5.231 -1.158  0.423  1.446  4.492 

Coefficients:
                             Estimate Std. Error t value Pr(>|t|)    
(Intercept)                  5.684779   0.396137  14.351  < 2e-16 ***
edad_primera_rsex           -0.017465   0.004425  -3.947 7.96e-05 ***
sexMujer                    -0.360791   0.030858 -11.692  < 2e-16 ***
orientacion_sexual_separado -0.052982   0.012477  -4.246 2.18e-05 ***
edad                        -0.026080   0.004952  -5.267 1.40e-07 ***
cohorteCohorte 1951-1960     0.502219   0.087169   5.761 8.49e-09 ***
cohorteCohorte 1961-1970     0.731319   0.120681   6.060 1.39e-09 ***
cohorteCohorte 1971-1980     0.861353   0.162963   5.286 1.27e-07 ***
cohorteCohorte 1981-1990     0.831484   0.210517   3.950 7.86e-05 ***
cohorteCohorte 1991-2000     0.597319   0.255378   2.339 0.019349 *  
cohorteCohorte 2001-2005     0.271129   0.294317   0.921 0.356951    
educSecundaria incompleta    0.170301   0.051936   3.279 0.001044 ** 
educSecondarian completa     0.388808   0.045655   8.516  < 2e-16 ***
educUniversitario            0.464984   0.047432   9.803  < 2e-16 ***
sexual_conversSi             0.091151   0.034890   2.613 0.008995 ** 
sexual_conversNo responde    0.032265   0.144041   0.224 0.822760    
sexual_educRegular           0.118845   0.035977   3.303 0.000957 ***
sexual_educBien              0.296057   0.039100   7.572 3.87e-14 ***
sexual_educSin respuesta     0.117389   0.067402   1.742 0.081592 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.856 on 16885 degrees of freedom
Multiple R-squared:  0.1374,    Adjusted R-squared:  0.1365 
F-statistic: 149.4 on 18 and 16885 DF,  p-value: < 2.2e-16

Click para ver el código

tidy(m1)

# A tibble: 19 × 5
   term                        estimate std.error statistic  p.value
   <chr>                          <dbl>     <dbl>     <dbl>    <dbl>
 1 (Intercept)                   5.68     0.396      14.4   1.98e-46
 2 edad_primera_rsex            -0.0175   0.00443    -3.95  7.96e- 5
 3 sexMujer                     -0.361    0.0309    -11.7   1.86e-31
 4 orientacion_sexual_separado  -0.0530   0.0125     -4.25  2.18e- 5
 5 edad                         -0.0261   0.00495    -5.27  1.40e- 7
 6 cohorteCohorte 1951-1960      0.502    0.0872      5.76  8.49e- 9
 7 cohorteCohorte 1961-1970      0.731    0.121       6.06  1.39e- 9
 8 cohorteCohorte 1971-1980      0.861    0.163       5.29  1.27e- 7
 9 cohorteCohorte 1981-1990      0.831    0.211       3.95  7.86e- 5
10 cohorteCohorte 1991-2000      0.597    0.255       2.34  1.93e- 2
11 cohorteCohorte 2001-2005      0.271    0.294       0.921 3.57e- 1
12 educSecundaria incompleta     0.170    0.0519      3.28  1.04e- 3
13 educSecondarian completa      0.389    0.0457      8.52  1.78e-17
14 educUniversitario             0.465    0.0474      9.80  1.25e-22
15 sexual_conversSi              0.0912   0.0349      2.61  8.99e- 3
16 sexual_conversNo responde     0.0323   0.144       0.224 8.23e- 1
17 sexual_educRegular            0.119    0.0360      3.30  9.57e- 4
18 sexual_educBien               0.296    0.0391      7.57  3.87e-14
19 sexual_educSin respuesta      0.117    0.0674      1.74  8.16e- 2

Click para ver el código

screenreg(l=list(m1), digits = 3)


==========================================
                             Model 1      
------------------------------------------
(Intercept)                      5.685 ***
                                (0.396)   
edad_primera_rsex               -0.017 ***
                                (0.004)   
sexMujer                        -0.361 ***
                                (0.031)   
orientacion_sexual_separado     -0.053 ***
                                (0.012)   
edad                            -0.026 ***
                                (0.005)   
cohorteCohorte 1951-1960         0.502 ***
                                (0.087)   
cohorteCohorte 1961-1970         0.731 ***
                                (0.121)   
cohorteCohorte 1971-1980         0.861 ***
                                (0.163)   
cohorteCohorte 1981-1990         0.831 ***
                                (0.211)   
cohorteCohorte 1991-2000         0.597 *  
                                (0.255)   
cohorteCohorte 2001-2005         0.271    
                                (0.294)   
educSecundaria incompleta        0.170 ** 
                                (0.052)   
educSecondarian completa         0.389 ***
                                (0.046)   
educUniversitario                0.465 ***
                                (0.047)   
sexual_conversSi                 0.091 ** 
                                (0.035)   
sexual_conversNo responde        0.032    
                                (0.144)   
sexual_educRegular               0.119 ***
                                (0.036)   
sexual_educBien                  0.296 ***
                                (0.039)   
sexual_educSin respuesta         0.117    
                                (0.067)   
------------------------------------------
R^2                              0.137    
Adj. R^2                         0.137    
Num. obs.                    16904        
==========================================
*** p < 0.001; ** p < 0.01; * p < 0.05

Click para ver el código

# Con diseño muestral
# Definir diseño muestral ENSSEX
enssex_design <- svydesign(
  id = ~varunit,   # UPM
  strata = ~varstrat,    # Estratos
  weights = ~ w_personas_cal,  # Ponderador muestral
  data = data,
  nest = TRUE,
)

options(survey.lonely.psu = "adjust")

# Estimamos
m2 <- svyglm(bienestar ~ edad_primera_rsex + sex + orientacion_sexual_separado + 
        edad + cohorte + educ + sexual_convers + sexual_educ, design = enssex_design)

screenreg(l=list(m1, m2), digits = 3)


=========================================================
                             Model 1        Model 2      
---------------------------------------------------------
(Intercept)                      5.685 ***      5.387 ***
                                (0.396)        (0.554)   
edad_primera_rsex               -0.017 ***     -0.030 ***
                                (0.004)        (0.007)   
sexMujer                        -0.361 ***     -0.409 ***
                                (0.031)        (0.048)   
orientacion_sexual_separado     -0.053 ***     -0.063 ** 
                                (0.012)        (0.019)   
edad                            -0.026 ***     -0.017 *  
                                (0.005)        (0.007)   
cohorteCohorte 1951-1960         0.502 ***      0.554 ***
                                (0.087)        (0.115)   
cohorteCohorte 1961-1970         0.731 ***      0.911 ***
                                (0.121)        (0.165)   
cohorteCohorte 1971-1980         0.861 ***      1.036 ***
                                (0.163)        (0.225)   
cohorteCohorte 1981-1990         0.831 ***      1.015 ***
                                (0.211)        (0.281)   
cohorteCohorte 1991-2000         0.597 *        0.814 *  
                                (0.255)        (0.344)   
cohorteCohorte 2001-2005         0.271          0.631    
                                (0.294)        (0.400)   
educSecundaria incompleta        0.170 **       0.268 ** 
                                (0.052)        (0.088)   
educSecondarian completa         0.389 ***      0.526 ***
                                (0.046)        (0.076)   
educUniversitario                0.465 ***      0.579 ***
                                (0.047)        (0.079)   
sexual_conversSi                 0.091 **       0.087 *  
                                (0.035)        (0.043)   
sexual_conversNo responde        0.032         -0.330    
                                (0.144)        (0.191)   
sexual_educRegular               0.119 ***      0.129 ** 
                                (0.036)        (0.049)   
sexual_educBien                  0.296 ***      0.279 ***
                                (0.039)        (0.058)   
sexual_educSin respuesta         0.117         -0.057    
                                (0.067)        (0.094)   
---------------------------------------------------------
R^2                              0.137                   
Adj. R^2                         0.137                   
Num. obs.                    16904          16904        
Deviance                                    51158.348    
Dispersion                                      3.027    
=========================================================
*** p < 0.001; ** p < 0.01; * p < 0.05

Regresión logística

Cuando la variable dependiente es binaria \(Y \in {0,1}\), usamos regresión logística ponderada:

Donde: - \(P(Y_i = 1)\) es la probabilidad de que ocurra el evento (Ej: uso de anticonceptivos = 1, no uso = 0). - \(\frac{P(Y_i = 1)}{1 - P(Y_i = 1)}\) es el odds ratio (razón de probabilidades). - \(\beta_p\) representa el efecto de cada variable independiente sobre la probabilidad del evento.

Click para ver el código

# Modelo sin diseño muestral
l1 <- glm(factor(bienestar>6) ~ edad_primera_rsex + sex + orientacion_sexual_separado + 
        edad + cohorte + educ + sexual_convers + sexual_educ,
        data=data,
        family = binomial(link="logit"))

# Modelo con diseño muestral
l2 <- svyglm(factor(bienestar>6) ~ edad_primera_rsex + sex + orientacion_sexual_separado + 
        edad + cohorte + educ + sexual_convers + sexual_educ,
       design = enssex_design, 
       family = quasibinomial())

# Reportamos
screenreg(l=list(l1, l2), digits = 3)


=========================================================
                             Model 1        Model 2      
---------------------------------------------------------
(Intercept)                     -0.595         -1.543 *  
                                (0.506)        (0.739)   
edad_primera_rsex               -0.019 **      -0.028 ** 
                                (0.006)        (0.009)   
sexMujer                        -0.017         -0.059    
                                (0.038)        (0.065)   
orientacion_sexual_separado     -0.030         -0.029    
                                (0.016)        (0.030)   
edad                            -0.016 *       -0.003    
                                (0.006)        (0.009)   
cohorteCohorte 1951-1960         0.286 *        0.448 *  
                                (0.133)        (0.197)   
cohorteCohorte 1961-1970         0.547 **       0.940 ***
                                (0.167)        (0.232)   
cohorteCohorte 1971-1980         0.711 ***      1.212 ***
                                (0.216)        (0.307)   
cohorteCohorte 1981-1990         0.720 **       1.229 ** 
                                (0.274)        (0.385)   
cohorteCohorte 1991-2000         0.641          1.236 ** 
                                (0.329)        (0.464)   
cohorteCohorte 2001-2005         0.509          1.248 *  
                                (0.377)        (0.535)   
educSecundaria incompleta        0.006          0.088    
                                (0.069)        (0.115)   
educSecondarian completa         0.063          0.213 *  
                                (0.059)        (0.099)   
educUniversitario                0.085          0.214 *  
                                (0.061)        (0.102)   
sexual_conversSi                 0.111 **       0.123    
                                (0.041)        (0.064)   
sexual_conversNo responde       -0.197         -0.534    
                                (0.179)        (0.296)   
sexual_educRegular              -0.031          0.029    
                                (0.044)        (0.065)   
sexual_educBien                  0.299 ***      0.384 ***
                                (0.046)        (0.079)   
sexual_educSin respuesta         0.160         -0.077    
                                (0.088)        (0.122)   
---------------------------------------------------------
AIC                          19346.677                   
BIC                          19493.647                   
Log Likelihood               -9654.338                   
Deviance                     19308.677      19842.472    
Num. obs.                    16904          16904        
Dispersion                                      1.001    
=========================================================
*** p < 0.001; ** p < 0.01; * p < 0.05

Click para ver el código

tidy(l1, exponentiate = TRUE)

# A tibble: 19 × 5
   term                        estimate std.error statistic  p.value
   <chr>                          <dbl>     <dbl>     <dbl>    <dbl>
 1 (Intercept)                    0.552   0.506     -1.17   2.40e- 1
 2 edad_primera_rsex              0.981   0.00597   -3.27   1.09e- 3
 3 sexMujer                       0.983   0.0376    -0.445  6.56e- 1
 4 orientacion_sexual_separado    0.970   0.0157    -1.93   5.32e- 2
 5 edad                           0.984   0.00628   -2.52   1.18e- 2
 6 cohorteCohorte 1951-1960       1.33    0.133      2.15   3.19e- 2
 7 cohorteCohorte 1961-1970       1.73    0.167      3.28   1.04e- 3
 8 cohorteCohorte 1971-1980       2.04    0.216      3.29   9.91e- 4
 9 cohorteCohorte 1981-1990       2.05    0.274      2.63   8.56e- 3
10 cohorteCohorte 1991-2000       1.90    0.329      1.95   5.14e- 2
11 cohorteCohorte 2001-2005       1.66    0.377      1.35   1.76e- 1
12 educSecundaria incompleta      1.01    0.0688     0.0827 9.34e- 1
13 educSecondarian completa       1.06    0.0594     1.06   2.90e- 1
14 educUniversitario              1.09    0.0610     1.39   1.63e- 1
15 sexual_conversSi               1.12    0.0409     2.71   6.67e- 3
16 sexual_conversNo responde      0.821   0.179     -1.10   2.71e- 1
17 sexual_educRegular             0.970   0.0441    -0.700  4.84e- 1
18 sexual_educBien                1.35    0.0460     6.51   7.67e-11
19 sexual_educSin respuesta       1.17    0.0881     1.82   6.91e- 2

Click para ver el código

tidy(l2, exponentiate = TRUE)

# A tibble: 19 × 5
   term                        estimate std.error statistic    p.value
   <chr>                          <dbl>     <dbl>     <dbl>      <dbl>
 1 (Intercept)                    0.214   0.739      -2.09  0.0377    
 2 edad_primera_rsex              0.973   0.00935    -2.94  0.00351   
 3 sexMujer                       0.942   0.0646     -0.917 0.360     
 4 orientacion_sexual_separado    0.971   0.0297     -0.981 0.328     
 5 edad                           0.997   0.00894    -0.315 0.753     
 6 cohorteCohorte 1951-1960       1.57    0.197       2.27  0.0238    
 7 cohorteCohorte 1961-1970       2.56    0.232       4.05  0.0000671 
 8 cohorteCohorte 1971-1980       3.36    0.307       3.95  0.0000990 
 9 cohorteCohorte 1981-1990       3.42    0.385       3.19  0.00158   
10 cohorteCohorte 1991-2000       3.44    0.464       2.67  0.00813   
11 cohorteCohorte 2001-2005       3.48    0.535       2.33  0.0205    
12 educSecundaria incompleta      1.09    0.115       0.765 0.445     
13 educSecondarian completa       1.24    0.0985      2.16  0.0313    
14 educUniversitario              1.24    0.102       2.09  0.0374    
15 sexual_conversSi               1.13    0.0637      1.93  0.0543    
16 sexual_conversNo responde      0.586   0.296      -1.80  0.0726    
17 sexual_educRegular             1.03    0.0653      0.450 0.653     
18 sexual_educBien                1.47    0.0788      4.87  0.00000188
19 sexual_educSin respuesta       0.926   0.122      -0.634 0.527

Aplicación I: análisis de sobrevivencia

El modelo de Cox (o regresión de riesgos proporcionales) se usa para tiempo hasta un evento (Ej: edad de primera relación sexual). Se expresa como:

\[h(t | X) = h_0(t) \exp(\beta_1 X_1 + \beta_2 X_2 + \dots + \beta_p X_p)\]

Donde: - \(h(t | X)\) es el riesgo instantáneo de que ocurra el evento en el tiempo ( t ). - \(h_0(t)\) es la función de riesgo base. - \(\beta_p\) son los coeficientes de las variables explicativas ( X_p ).

Click para ver el código

# Modelos de COX
cox <- coxph(Surv(edad_primera_rsex) ~ bienestar + sex + orientacion_sexual_separado + edad + cohorte + educ + sexual_convers + sexual_educ, data = data)

cox_ponderado <- svycoxph(Surv(edad_primera_rsex) ~ bienestar + sex + orientacion_sexual_separado + 
        edad + cohorte + educ + sexual_convers + sexual_educ, design = enssex_design)

# Modelos de AFT
aft1 <- survreg(Surv(edad_primera_rsex) ~bienestar + sex + orientacion_sexual_separado + 
        edad + cohorte + educ + sexual_convers + sexual_educ, data = data, dist = "gaussian")
aft2 <- survreg(Surv(edad_primera_rsex) ~bienestar + sex + orientacion_sexual_separado + 
        edad + cohorte + educ + sexual_convers + sexual_educ, data = data, dist = "logistic")
aft3 <- survreg(Surv(edad_primera_rsex) ~bienestar + sex + orientacion_sexual_separado + 
        edad + cohorte + educ + sexual_convers + sexual_educ, data = data, dist = "extreme")
aft4 <- survreg(Surv(edad_primera_rsex) ~bienestar + sex + orientacion_sexual_separado + 
        edad + cohorte + educ + sexual_convers + sexual_educ, data = data, dist = "weibull")
aft5 <- survreg(Surv(edad_primera_rsex) ~bienestar + sex + orientacion_sexual_separado + 
        edad + cohorte + educ + sexual_convers + sexual_educ, data = data, dist = "exponential")
aft6 <- survreg(Surv(edad_primera_rsex) ~bienestar + sex + orientacion_sexual_separado + 
        edad + cohorte + educ + sexual_convers + sexual_educ, data = data, dist = "lognormal")

# Reportamos 
summary(cox)

Call:
coxph(formula = Surv(edad_primera_rsex) ~ bienestar + sex + orientacion_sexual_separado + 
    edad + cohorte + educ + sexual_convers + sexual_educ, data = data)

  n= 16904, number of events= 16904 

                                 coef exp(coef)  se(coef)       z Pr(>|z|)    
bienestar                    0.017808  1.017967  0.004111   4.332 1.48e-05 ***
sexMujer                    -0.388540  0.678046  0.016454 -23.614  < 2e-16 ***
orientacion_sexual_separado  0.011343  1.011408  0.006704   1.692  0.09065 .  
edad                        -0.024532  0.975766  0.002706  -9.065  < 2e-16 ***
cohorteCohorte 1951-1960    -0.134674  0.874001  0.046993  -2.866  0.00416 ** 
cohorteCohorte 1961-1970    -0.314725  0.729990  0.065426  -4.810 1.51e-06 ***
cohorteCohorte 1971-1980    -0.260730  0.770489  0.088389  -2.950  0.00318 ** 
cohorteCohorte 1981-1990    -0.260043  0.771018  0.114083  -2.279  0.02264 *  
cohorteCohorte 1991-2000    -0.252746  0.776665  0.138491  -1.825  0.06800 .  
cohorteCohorte 2001-2005    -0.107290  0.898265  0.159751  -0.672  0.50183    
educSecundaria incompleta   -0.077818  0.925133  0.028045  -2.775  0.00552 ** 
educSecondarian completa    -0.180989  0.834445  0.024533  -7.377 1.61e-13 ***
educUniversitario           -0.314788  0.729944  0.025438 -12.375  < 2e-16 ***
sexual_conversSi            -0.015008  0.985104  0.018770  -0.800  0.42395    
sexual_conversNo responde    0.007278  1.007305  0.077827   0.094  0.92549    
sexual_educRegular          -0.038835  0.961910  0.019435  -1.998  0.04570 *  
sexual_educBien             -0.059724  0.942024  0.021169  -2.821  0.00478 ** 
sexual_educSin respuesta    -0.028206  0.972188  0.036386  -0.775  0.43823    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

                            exp(coef) exp(-coef) lower .95 upper .95
bienestar                      1.0180     0.9823    1.0098    1.0262
sexMujer                       0.6780     1.4748    0.6565    0.7003
orientacion_sexual_separado    1.0114     0.9887    0.9982    1.0248
edad                           0.9758     1.0248    0.9706    0.9810
cohorteCohorte 1951-1960       0.8740     1.1442    0.7971    0.9583
cohorteCohorte 1961-1970       0.7300     1.3699    0.6421    0.8299
cohorteCohorte 1971-1980       0.7705     1.2979    0.6479    0.9162
cohorteCohorte 1981-1990       0.7710     1.2970    0.6165    0.9642
cohorteCohorte 1991-2000       0.7767     1.2876    0.5920    1.0189
cohorteCohorte 2001-2005       0.8983     1.1133    0.6568    1.2285
educSecundaria incompleta      0.9251     1.0809    0.8757    0.9774
educSecondarian completa       0.8344     1.1984    0.7953    0.8755
educUniversitario              0.7299     1.3700    0.6944    0.7673
sexual_conversSi               0.9851     1.0151    0.9495    1.0220
sexual_conversNo responde      1.0073     0.9927    0.8648    1.1733
sexual_educRegular             0.9619     1.0396    0.9260    0.9993
sexual_educBien                0.9420     1.0615    0.9037    0.9819
sexual_educSin respuesta       0.9722     1.0286    0.9053    1.0441

Concordance= 0.647  (se = 0.003 )
Likelihood ratio test= 2755  on 18 df,   p=<2e-16
Wald test            = 2838  on 18 df,   p=<2e-16
Score (logrank) test = 2939  on 18 df,   p=<2e-16

Click para ver el código

summary(cox_ponderado)

Stratified 1 - level Cluster Sampling design (with replacement)
With (372) clusters.
svydesign(id = ~varunit, strata = ~varstrat, weights = ~w_personas_cal, 
    data = data, nest = TRUE, )

Call:
svycoxph(formula = Surv(edad_primera_rsex) ~ bienestar + sex + 
    orientacion_sexual_separado + edad + cohorte + educ + sexual_convers + 
    sexual_educ, design = enssex_design)

  n= 16904, number of events= 16904 

                                  coef  exp(coef)   se(coef)  robust se       z
bienestar                    3.159e-02  1.032e+00  4.331e-03  6.418e-03   4.923
sexMujer                    -3.512e-01  7.039e-01  1.571e-02  2.544e-02 -13.803
orientacion_sexual_separado -4.932e-03  9.951e-01  8.071e-03  1.824e-02  -0.270
edad                        -2.475e-02  9.756e-01  2.696e-03  4.252e-03  -5.819
cohorteCohorte 1951-1960    -1.444e-01  8.655e-01  4.943e-02  7.105e-02  -2.033
cohorteCohorte 1961-1970    -3.588e-01  6.985e-01  6.747e-02  9.497e-02  -3.778
cohorteCohorte 1971-1980    -2.749e-01  7.597e-01  9.050e-02  1.379e-01  -1.994
cohorteCohorte 1981-1990    -3.729e-01  6.887e-01  1.152e-01  1.830e-01  -2.038
cohorteCohorte 1991-2000    -4.382e-01  6.452e-01  1.397e-01  2.146e-01  -2.041
cohorteCohorte 2001-2005    -2.968e-01  7.432e-01  1.612e-01  2.522e-01  -1.177
educSecundaria incompleta   -9.524e-02  9.092e-01  3.231e-02  4.782e-02  -1.992
educSecondarian completa    -2.573e-01  7.731e-01  2.617e-02  5.084e-02  -5.062
educUniversitario           -3.366e-01  7.142e-01  2.614e-02  4.758e-02  -7.074
sexual_conversSi            -5.884e-05  9.999e-01  1.852e-02  2.960e-02  -0.002
sexual_conversNo responde   -2.425e-02  9.760e-01  8.717e-02  1.090e-01  -0.223
sexual_educRegular          -2.107e-02  9.791e-01  1.951e-02  3.492e-02  -0.603
sexual_educBien              1.586e-02  1.016e+00  2.102e-02  3.351e-02   0.473
sexual_educSin respuesta    -7.453e-02  9.282e-01  3.196e-02  5.383e-02  -1.384
                            Pr(>|z|)    
bienestar                   8.52e-07 ***
sexMujer                     < 2e-16 ***
orientacion_sexual_separado 0.786901    
edad                        5.91e-09 ***
cohorteCohorte 1951-1960    0.042092 *  
cohorteCohorte 1961-1970    0.000158 ***
cohorteCohorte 1971-1980    0.046185 *  
cohorteCohorte 1981-1990    0.041569 *  
cohorteCohorte 1991-2000    0.041202 *  
cohorteCohorte 2001-2005    0.239312    
educSecundaria incompleta   0.046404 *  
educSecondarian completa    4.15e-07 ***
educUniversitario           1.51e-12 ***
sexual_conversSi            0.998414    
sexual_conversNo responde   0.823862    
sexual_educRegular          0.546203    
sexual_educBien             0.636021    
sexual_educSin respuesta    0.166217    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

                            exp(coef) exp(-coef) lower .95 upper .95
bienestar                      1.0321     0.9689    1.0192    1.0452
sexMujer                       0.7039     1.4207    0.6696    0.7399
orientacion_sexual_separado    0.9951     1.0049    0.9601    1.0313
edad                           0.9756     1.0251    0.9675    0.9837
cohorteCohorte 1951-1960       0.8655     1.1554    0.7530    0.9949
cohorteCohorte 1961-1970       0.6985     1.4316    0.5799    0.8414
cohorteCohorte 1971-1980       0.7597     1.3164    0.5798    0.9954
cohorteCohorte 1981-1990       0.6887     1.4520    0.4811    0.9859
cohorteCohorte 1991-2000       0.6452     1.5499    0.4236    0.9827
cohorteCohorte 2001-2005       0.7432     1.3455    0.4533    1.2184
educSecundaria incompleta      0.9092     1.0999    0.8278    0.9985
educSecondarian completa       0.7731     1.2935    0.6998    0.8541
educUniversitario              0.7142     1.4001    0.6506    0.7840
sexual_conversSi               0.9999     1.0001    0.9436    1.0597
sexual_conversNo responde      0.9760     1.0245    0.7884    1.2084
sexual_educRegular             0.9791     1.0213    0.9144    1.0485
sexual_educBien                1.0160     0.9843    0.9514    1.0850
sexual_educSin respuesta       0.9282     1.0774    0.8352    1.0315

Concordance= 0.629  (se = 0.004 )
Likelihood ratio test= NA  on 18 df,   p=NA
Wald test            = 1197  on 18 df,   p=<2e-16
Score (logrank) test = NA  on 18 df,   p=NA

  (Note: the likelihood ratio and score tests assume independence of
     observations within a cluster, the Wald and robust score tests do not).

Click para ver el código

screenreg(l=list(cox, aft1, aft2, aft3, aft4, aft5, aft6), digits = 3,
        custom.model.names = c("Cox", "AFT Gaussiano", "AFT Logístico", "AFT Extremo", 
                               "AFT Weibull", "AFT Exponential", "AFT lognormal"))


============================================================================================================================================
                             Cox             AFT Gaussiano   AFT Logístico   AFT Extremo     AFT Weibull     AFT Exponential  AFT lognormal 
--------------------------------------------------------------------------------------------------------------------------------------------
bienestar                         0.018 ***      -0.053 ***      -0.050 ***      -0.082 ***      -0.006 ***      -0.003           -0.003 ***
                                 (0.004)         (0.013)         (0.011)         (0.022)         (0.001)         (0.004)          (0.001)   
sexMujer                         -0.389 ***       1.449 ***       1.264 ***       2.076 ***       0.087 ***       0.080 ***        0.080 ***
                                 (0.016)         (0.053)         (0.044)         (0.088)         (0.003)         (0.016)          (0.003)   
orientacion_sexual_separado       0.011          -0.045 *        -0.033          -0.046          -0.003 *        -0.002           -0.002 *  
                                 (0.007)         (0.022)         (0.018)         (0.037)         (0.001)         (0.007)          (0.001)   
edad                             -0.025 ***       0.077 ***       0.063 ***       0.250 ***       0.007 ***       0.004            0.004 ***
                                 (0.003)         (0.009)         (0.007)         (0.013)         (0.001)         (0.003)          (0.000)   
cohorteCohorte 1951-1960         -0.135 **        0.305 *         0.258           2.100 ***       0.008           0.018            0.019 *  
                                 (0.047)         (0.152)         (0.133)         (0.258)         (0.010)         (0.047)          (0.008)   
cohorteCohorte 1961-1970         -0.315 ***       0.738 ***       0.630 ***       3.835 ***       0.042 **        0.043            0.042 ***
                                 (0.065)         (0.210)         (0.181)         (0.346)         (0.014)         (0.065)          (0.011)   
cohorteCohorte 1971-1980         -0.261 **        0.432           0.284           4.483 ***       0.042 *         0.028            0.026    
                                 (0.088)         (0.283)         (0.242)         (0.460)         (0.019)         (0.088)          (0.014)   
cohorteCohorte 1981-1990         -0.260 *         0.451           0.301           5.053 ***       0.042           0.028            0.026    
                                 (0.114)         (0.366)         (0.312)         (0.591)         (0.024)         (0.114)          (0.018)   
cohorteCohorte 1991-2000         -0.253           0.574           0.366           6.534 ***       0.060 *         0.033            0.029    
                                 (0.138)         (0.444)         (0.378)         (0.712)         (0.029)         (0.138)          (0.022)   
cohorteCohorte 2001-2005         -0.107           0.655           0.449           7.526 ***       0.061           0.034            0.031    
                                 (0.160)         (0.511)         (0.434)         (0.816)         (0.034)         (0.159)          (0.026)   
educSecundaria incompleta        -0.078 **        0.437 ***       0.477 ***       0.364 *         0.010           0.023            0.024 ***
                                 (0.028)         (0.090)         (0.078)         (0.151)         (0.006)         (0.028)          (0.005)   
educSecondarian completa         -0.181 ***       0.844 ***       0.864 ***       0.873 ***       0.023 ***       0.045            0.047 ***
                                 (0.025)         (0.079)         (0.069)         (0.131)         (0.005)         (0.025)          (0.004)   
educUniversitario                -0.315 ***       1.259 ***       1.271 ***       0.949 ***       0.043 ***       0.068 **         0.071 ***
                                 (0.025)         (0.082)         (0.071)         (0.136)         (0.005)         (0.025)          (0.004)   
sexual_conversSi                 -0.015           0.067           0.012           0.576 ***       0.007           0.003            0.003    
                                 (0.019)         (0.061)         (0.050)         (0.099)         (0.004)         (0.019)          (0.003)   
sexual_conversNo responde         0.007          -0.000           0.095          -0.467          -0.018           0.001            0.003    
                                 (0.078)         (0.250)         (0.205)         (0.412)         (0.016)         (0.078)          (0.013)   
sexual_educRegular               -0.039 *         0.162 **        0.146 **        0.609 ***       0.020 ***       0.009            0.008 ** 
                                 (0.019)         (0.063)         (0.053)         (0.103)         (0.004)         (0.019)          (0.003)   
sexual_educBien                  -0.060 **        0.182 **        0.187 **        0.274 *         0.031 ***       0.010            0.010 ** 
                                 (0.021)         (0.068)         (0.057)         (0.113)         (0.004)         (0.021)          (0.003)   
sexual_educSin respuesta         -0.028           0.115          -0.031           1.514 ***       0.056 ***       0.006            0.002    
                                 (0.036)         (0.117)         (0.102)         (0.196)         (0.008)         (0.036)          (0.006)   
(Intercept)                                      12.279 ***      12.863 ***       1.677           2.540 ***       2.573 ***        2.574 ***
                                                 (0.686)         (0.584)         (1.074)         (0.045)         (0.213)          (0.035)   
Log(scale)                                        1.171 ***       0.465 ***       1.666 ***      -1.560 ***                       -1.816 ***
                                                 (0.005)         (0.007)         (0.004)         (0.005)                          (0.005)   
--------------------------------------------------------------------------------------------------------------------------------------------
AIC                          292616.244       87594.104       84196.554      102673.920       91600.274      130979.129        83268.230    
R^2                               0.150                                                                                                     
Max. R^2                          1.000                                                                                                     
Num. events                   16904                                                                                                         
Num. obs.                     16904           16904           16904           16904           16904           16904            16904        
Missings                          0                                                                                                         
PH test                           0.000                                                                                                     
BIC                                           87748.810       84351.260      102828.626       91754.980      131126.100        83422.936    
Log Likelihood                               -43777.052      -42078.277      -51316.960      -45780.137      -65470.564       -41614.115    
============================================================================================================================================
*** p < 0.001; ** p < 0.01; * p < 0.05

Aplicación II: Análisis de clusters y segmentación

Agrupa observaciones según similitudes en sus características:

\[\arg\min_S \sum_{i=1}^{k} \sum_{x \in S_i} | x - \mu_i |^2\]

Donde: - \(S_i\) son los grupos (clusters). - \(\mu_i\) es el centroide de cada grupo. - \(x\) son los valores de las observaciones.

Para esto utilizaremos un set de preguntas:

Click para ver el código

# Usemos un set de preguntas
data <- rio::import("https://datos.gob.cl/dataset/c6983439-49f6-4e71-85fe-e8de6e73dae0/resource/ed81f50c-1c7d-43d9-9083-dfc161e0cd66/download/20240516_enssex_data.rdata") |> 
  select(ends_with("p9")) |> 
  mutate(across(everything(), ~ zap_labels(na_if(., 9)))) |> 
  drop_na()

Click para ver el código

glimpse(data)

Rows: 18,995
Columns: 6
$ i_1_p9 <dbl> 7, 4, 4, 6, 1, 4, 3, 7, 7, 1, 1, 3, 1, 3, 7, 7, 4, 7, 5, 5, 4, …
$ i_2_p9 <dbl> 6, 4, 4, 4, 1, 5, 7, 7, 7, 4, 4, 4, 3, 3, 5, 7, 6, 7, 6, 7, 3, …
$ i_3_p9 <dbl> 7, 3, 5, 5, 1, 6, 5, 7, 5, 2, 7, 2, 4, 4, 6, 7, 6, 7, 4, 5, 6, …
$ i_4_p9 <dbl> 6, 6, 4, 6, 3, 6, 5, 6, 7, 6, 1, 1, 7, 4, 5, 5, 3, 7, 5, 5, 6, …
$ i_5_p9 <dbl> 7, 6, 5, 4, 2, 5, 5, 6, 7, 7, 3, 6, 6, 5, 6, 7, 5, 6, 5, 5, 7, …
$ i_6_p9 <dbl> 7, 6, 5, 4, 6, 4, 5, 6, 5, 1, 4, 2, 6, 5, 6, 7, 6, 7, 4, 5, 4, …

Veamos los descriptivos (sin ponderar):

Click para ver el código

st(data,
   digits = 3, 
   out="kable", 
   add.median = TRUE,
   fixed.digits = TRUE, 
   simple.kable = FALSE,
   title="",
   numformat = NA) %>%
  kableExtra::kable_styling(latex_options = c("scale_down", "hold_position"),
                            full_width = FALSE, fixed_thead = T)

Variable	N	Mean	Std. Dev.	Min	Pctl. 25	Pctl. 50	Pctl. 75	Max
i_1_p9	18995	5.761	1.363	1	5	6	7	7
i_2_p9	18995	5.513	1.326	1	5	6	7	7
i_3_p9	18995	5.135	1.922	1	4	6	7	7
i_4_p9	18995	5.213	1.554	1	4	5	6	7
i_5_p9	18995	5.926	1.284	1	5	6	7	7
i_6_p9	18995	4.855	2.057	1	4	5	7	7

¿Cuántos cluster?

Click para ver el código

# Importante: estandarizar
data_zscore <- scale(data)

fviz_nbclust(data_zscore, 
  kmeans, 
  k.max = 6,
  method = "wss") +  
  theme_minimal()

Aplicamos nuestra estrategia de clustering

Click para ver el código

set.seed(123)
model.km <- kmeans(data_zscore, 
                   iter.max=5,
                   centers=3, 
                   nstart=100, 
                   trace = FALSE)

summary(model.km)

             Length Class  Mode   
cluster      18995  -none- numeric
centers         18  -none- numeric
totss            1  -none- numeric
withinss         3  -none- numeric
tot.withinss     1  -none- numeric
betweenss        1  -none- numeric
size             3  -none- numeric
iter             1  -none- numeric
ifault           1  -none- numeric

Click para ver el código

# Veamos los cluster
data$cluster_kmedias <- model.km$cluster
head(data)

  i_1_p9 i_2_p9 i_3_p9 i_4_p9 i_5_p9 i_6_p9 cluster_kmedias
1      7      6      7      6      7      7               3
2      4      4      3      6      6      6               2
3      4      4      5      4      5      5               2
4      6      4      5      6      4      4               2
5      1      1      1      3      2      6               2
6      4      5      6      6      5      4               2

Exploremos las categorías obtenidas:

Click para ver el código

# Visualicemos los cluster 
fviz_cluster(model.km, data = data, geom = "point") +
  labs(title=NULL) +
  theme_light()

Click para ver el código

# Exploremos los cluster
aggregate(data[1:6], by=list(model.km$cluster), mean)

  Group.1   i_1_p9   i_2_p9   i_3_p9   i_4_p9   i_5_p9   i_6_p9
1       1 6.161019 5.594303 2.696252 4.794903 6.128936 2.057871
2       2 4.331965 4.116706 4.154528 3.824941 4.523882 3.974065
3       3 6.239917 6.073779 6.284148 5.921961 6.452280 6.070788

Ejercicio propuesto

Piense en una pregunta de investigación que pueda ser respondida con datos con de la encuesta ENSSEX y proponga de forma intuitiva un modelo para este fenómeno. Implemente este modelo con las herramientas vistas en el workshop.